Introduction

This exercise was modified based on https://geodacenter.github.io/aot-workshop/Part2-AOT.html.

Goals for this lab: * Deepen the understanding of surface interpolation and the procedure of ordinary kriging. * Practice using R’s variogram and kriging functions (more at: https://gsp.humboldt.edu/olm/R/04_01_Variograms.html) to perform ordinary kriging for temperature point data.

Go through the steps below and understand the meaning of each script line. Then, answer the question at the end of this file.

1 Environment setup

library(lubridate) #data wrangling
library(sf)
library(sp) #spatial data wrangling & analysis
library(raster) #spatial raster data wrangling
library(terra) #spatial raster data wrangling
#install.packages("gstat")
library(gstat) #kriging and geostatistics
library(tmap) #modern data visualizations
library(leaflet) #modern data visualizations

# Define your working directory
setwd("/Users/jonathanbernard/Desktop/UIUC/Spring 2024/GGIS 224/GGIS224/Labs/Lab11")

2 Import temperature point data

nodes <- read.csv("Data_Lab11/AOT_node_temps.txt") #Read the file as csv
head(nodes) #Inspect the data: 31 points with temperature values in both Celsius and Fahrenheit degrees.
##        node_id avg_temp avg_temp_f  project_id  vsn
## 1 001e06109416 26.01265   78.82278 AoT_Chicago 090B
## 2 001e0610b9e9 25.87948   78.58306 AoT_Chicago   80
## 3 001e0610ba8f 28.30287   82.94517 AoT_Chicago  00D
## 4 001e0610bbf9 25.96714   78.74086 AoT_Chicago   20
## 5 001e0610bbff 27.93934   82.29081 AoT_Chicago   25
## 6 001e0610bc10 25.51712   77.93082 AoT_Chicago  01F
##                                  address      lat       lon         description
## 1         Elston and Cortland Chicago IL 41.91659 -87.66641 AoT Chicago (S) [C]
## 2 Broadway Ave & Lawrence Ave Chicago IL 41.96904 -87.65967 AoT Chicago (S) [C]
## 3           Cornell & 47th St Chicago IL 41.81034 -87.59023     AoT Chicago (S)
## 4       Western Ave & 69th St Chicago IL 41.76832 -87.68340 AoT Chicago (S) [C]
## 5       Western Ave & 18th St Chicago IL 41.85780 -87.68581     AoT Chicago (S)
## 6             State St & 87th Chicago IL 41.73631 -87.62418 AoT Chicago (S) [C]
##   start_timestamp
## 1     1/1/18 0:00
## 2    2/14/18 0:00
## 3     8/8/17 0:00
## 4    2/13/18 0:00
## 5   12/15/17 0:00
## 6    2/22/18 0:00

3 Convert to spatial data format

Convert the completed node data to spatial object format for plotting and spatial analytics.

node.temps <- nodes #Duplicate 'nodes' to another data file
coordinates(node.temps) <- node.temps[,c("lon", "lat")]
proj4string(node.temps) <- CRS("+init=epsg:4326")
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.

4 Inspect the converted data

There are 31 sensor points with temperature measurements in both Celsius and Fahrenheit degrees.

length(node.temps)
## [1] 31
head(node.temps) #avg_temp for Celsius; avg_temp_f for Fahrenheit. 
##        node_id avg_temp avg_temp_f  project_id  vsn
## 1 001e06109416 26.01265   78.82278 AoT_Chicago 090B
## 2 001e0610b9e9 25.87948   78.58306 AoT_Chicago   80
## 3 001e0610ba8f 28.30287   82.94517 AoT_Chicago  00D
## 4 001e0610bbf9 25.96714   78.74086 AoT_Chicago   20
## 5 001e0610bbff 27.93934   82.29081 AoT_Chicago   25
## 6 001e0610bc10 25.51712   77.93082 AoT_Chicago  01F
##                                  address      lat       lon         description
## 1         Elston and Cortland Chicago IL 41.91659 -87.66641 AoT Chicago (S) [C]
## 2 Broadway Ave & Lawrence Ave Chicago IL 41.96904 -87.65967 AoT Chicago (S) [C]
## 3           Cornell & 47th St Chicago IL 41.81034 -87.59023     AoT Chicago (S)
## 4       Western Ave & 69th St Chicago IL 41.76832 -87.68340 AoT Chicago (S) [C]
## 5       Western Ave & 18th St Chicago IL 41.85780 -87.68581     AoT Chicago (S)
## 6             State St & 87th Chicago IL 41.73631 -87.62418 AoT Chicago (S) [C]
##   start_timestamp
## 1     1/1/18 0:00
## 2    2/14/18 0:00
## 3     8/8/17 0:00
## 4    2/13/18 0:00
## 5   12/15/17 0:00
## 6    2/22/18 0:00

5 Plot temperature data

Confirm the success of spatial object transformation by simple plotting.

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(node.temps) + tm_dots()

6 Overlay community areas

We here add the Chicago Community Area spatial data to provide additional context for our maps. The Chicago Community Area shapefile is downloaded from the Chicago Data Portal website: https://data.cityofchicago.org/Facilities-Geographic-Boundaries/Boundaries-Community-Areas-current-/cauq-8yn6.

chiCA <- read_sf('Data_Lab11/ChicagoCommunityAreas/geo_export_fcf7548c-1d07-4c42-99a3-156e428257cc.shp')

Plot the temperature by sensor points, adding Community Areas as a background.

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chiCA) + tm_borders() + 
  tm_shape(node.temps) + tm_dots(col="avg_temp_f",size=0.3,title="average temp (F)")

7 Interpolate a temperature surface

We will use (semi)variograms to model the distribution of temperature data. We will then generate a grid on top of the Chicago area, and with the appropriate variogram model selected, use kriging to predict a temperature surface.

A variogram is a function that describes the degree of spatial autocorrelation across data. The final model uses the measure of variability between points at various distances. Points nearby are likely to have more similar values, and as distance between points increases, there is less likely to be similar values between points. In this application, we assume that temperature measurements that are further apart will vary more than measurements taken close together (so, a positive autocorrelation).

The variogram clearly has an outlier, though it may not influence our final predicted surface because the variogram will be fitted by a model.

tmp.vgm <- variogram(node.temps$avg_temp_f ~ 1, node.temps) #here variable avg_temp_f (temperature measurements in Fahrenheit) is used. 
plot(tmp.vgm)

8 Semivariogram fitting: using the Spherical model

We will generate a theoretical spherical model to approximate (fit) the experimental semivariogram.

tmp.fit.sph<- fit.variogram(tmp.vgm, model=vgm("Sph"))
plot(tmp.vgm, tmp.fit.sph)

9 Generate grid

Next, we’ll create a grid from the Chicago area. The following function will generate a n-by-n grid from a provided spatial data frame.

pt2grid <- function(ptframe,n) {
  bb <- bbox(ptframe)  
  ptcrs <- proj4string(ptframe)  
  xrange <- abs(bb[1,1] - bb[1,2])  
  yrange <- abs(bb[2,1] - bb[2,2])  
  cs <- c(xrange/n,yrange/n)  
  cc <- bb[,1] + (cs/2)  
  dc <- c(n,n)  
  x1 <- GridTopology(cellcentre.offset=cc,cellsize=cs,cells.dim=dc)  
  x2 <- SpatialGrid(grid=x1,proj4string=CRS(ptcrs))
  return(x2)
}

Let’s first generate a grid of 30 by 30 cells using the Chicago Community Area extent. Plot the grid for exploration.

chiCA = as_Spatial(chiCA) #converting chiCA (sf) to a spatial data frame object (needed for executing the pt2grid function)
chi.grid <- pt2grid(chiCA,30)
plot(chi.grid)

To get an even finer resolution, we generate a finer resolution of grid of 100 by 100 cells. Let’s use the finer grid for the following analysis.

chi.grid <- pt2grid(chiCA,100)
plot(chi.grid)

10 Prepare data for kriging

First, we make sure that all our data is in the same projection.

projection(chi.grid) <- crs("+init=epsg:4326")  
projection(node.temps) <-  crs("+init=epsg:4326")
projection(chiCA) <- crs("+init=epsg:4326")

Krige the data using the spherical model.

temp.kriged <- krige(node.temps$avg_temp_f ~ 1, node.temps, chi.grid, model = tmp.fit.sph) #Again, variable avg_temp_f is used for kriging. 
## [using ordinary kriging]
plot(temp.kriged)

Convert the files to SpatRaster and SpatVector files, clip to Chicago boundaries, and plot to confirm everything looks good.

chiCA.spat <- terra::vect(chiCA) #convert so we can use in Terra
temp.kriged.spat <- terra::rast(temp.kriged) #convert so we can use in Terra
chi.temp.kriged <- terra::crop(temp.kriged.spat, chiCA.spat, mask=TRUE) #use chiCA as a mask so we only see temperatures within the Chicago boundary
plot(chi.temp.kriged$var1.pred)
plot(node.temps, add = TRUE)
plot(chiCA, add = TRUE)

Map the kriged Chicago-area surface.

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chiCA) + tm_borders() + 
  tm_shape(node.temps) + tm_dots(size=0.01) +
  tm_shape(chi.temp.kriged$var1.pred) + tm_raster("var1.pred", style = "jenks", title = "Temperature (F)", palette = "BuPu") +
  tm_legend(position = c("left", "bottom")) #var1.pred is the ordinary kriging predictions, and #var1.var is the ordinary kriging variance
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.

Question

Now it is your turn. Try to krige the temperature values (this time, in Celsius degree) using a exponential fitting model, and map both of the predicted surface and its variance. Hint: to fit the semivariogram by an exponential model, change the argument of model in fit.variogram to “Exp”. Also see more: https://search.r-project.org/CRAN/refmans/gstat/html/fit.variogram.html.

# semivariogram fitting using the exponential model
tmp.vgm <- variogram(node.temps$avg_temp ~ 1, node.temps) 
tmp.fit.exp <- fit.variogram(tmp.vgm, model=vgm("Exp"))
plot(tmp.vgm, tmp.fit.exp)

# generating the grid
chi.grid <- pt2grid(chiCA,100)

# preparing the data for kriging analysis
projection(chi.grid) <- crs("+init=epsg:4326")  
projection(node.temps) <-  crs("+init=epsg:4326")
projection(chiCA) <- crs("+init=epsg:4326")

# krige the data using the exponential model
temp.kriged <- krige(node.temps$avg_temp ~ 1, node.temps, chi.grid, model = tmp.fit.exp)
## [using ordinary kriging]
# Converting the data to to SpatRaster and SpatVector
chiCA.spat <- terra::vect(chiCA)
temp.kriged.spat <- terra::rast(temp.kriged)
chi.temp.kriged <- terra::crop(temp.kriged.spat, chiCA.spat, mask=TRUE)

# mapping the kriged predicted surface using the view mode
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(chiCA) + tm_borders() + 
  tm_shape(node.temps) + tm_dots(size=0.01) +
  tm_shape(chi.temp.kriged$var1.pred) + tm_raster("var1.pred", style = "jenks", 
            title = "Temperature (C)", palette = "BuPu") +
  tm_legend(position = c("left", "bottom"))
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
## Use "fisher" instead of "jenks" for larger data sets
# mapping the kriging variance also using the view mode
tmap_mode("view") 
## tmap mode set to interactive viewing
tm_shape(chiCA) + tm_borders() +
  tm_shape(node.temps) + tm_dots(size=0.01) + 
  tm_shape(chi.temp.kriged$var1.var) + tm_raster("var1.var", style = "jenks",
            title = "Kriging Variance", palette = "Reds") +
  tm_legend(position = c("left", "bottom"))
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
## Use "fisher" instead of "jenks" for larger data sets

.